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We determine the zero temperature quantum phase diagram of a p x + ip y pairing model based 
on the exactly solvable hyperbolic Richardson-Gaudin model. We present analytical and large-scale 
numerical results for this model. In the continuum limit, the exact solution exhibits a third-order 
quantum phase transition, separating a strong-pairing from a weak-pairing phase. The mean held 
solution allows to connect these results to other models with p x + ip y pairing order. We define 
an experimentally accessible characteristic length scale, associated with the size of the Cooper 
pairs, that diverges at the transition point, indicating that the phase transition is of a confinement- 
deconfinement type without local order parameter. We show that this phase transition is not limited 
to the p x + ip y pairing model, but can be found in any representation of the hyperbolic Richardson- 
Gaudin model and is related to a symmetry that is absent in the rational Richardson-Gaudin model. 

PACS numbers: 74.20. Fg, 74.20.Rp, 71.10.Hf, 73.43.Nq 

I. INTRODUCTION 

One of the striking features of degenerate Fermi gases whose constituents interact through an attractive potential 
is that they can exhibit superfluidity or superconductivity. The constituents may be neutral atoms, as is the case in 
trapped Fermi gases of 40 K and 6 LP or superfiuid 3 He^, or charged electrons as in conventional metallic supercon- 
ductors, nucleons in heavy nuclei or nuclear matter^. The particular nature and symmetry of the attractive pairing 
interaction may lead to exotic superfiuid phases with complex order parameters, as is well-known in the case of liq- 
uid 3 He. For order parameters with certain rotational symmetries, a change in the coupling strength may induce a 
quantum phase transition (QPT) separating different kinds of superfiuid phases. In a Fermi system whose pairing 
interaction has an s-wave character, it is well understood that by increasing the coupling strength there is a crossover, 
and not a QPT, between a weak coupling Bardeen-Cooper-Schrieffer 4 (BCS) and a Bose-Einstein condensate (BEC) 
phase 5 . However, for higher rotational order, the ground state can exhibit a QPT between qualitatively different 
superfiuid states, with a corresponding non-analyticity in the ground state energy. In the present paper we focus on 
p-wave pairing that might be encountered in ultracold Fermi gases^Hor in exotic superconductors such as Sr 2 RuOj^. 

As a schematic model for such systems we study the quantum phase diagram of a spinless Fermi gas with p x + ip y 
pairing interaction symmetry 

k 2 

H P X +iPy = Y C fe C fc ~ V kk'clcl k C_ k ,C k ,, (1) 

k k,k' 

where Vkk' — 9k9k' 1S a separable interaction, with gk representing a complex function of wavevector k and symmetry 
Px + ipy Though schematic, this model captures the basic physics of the BEC and BCS sides of the phase diagram 
and offers detailed insights into the phase transition mechanism. Previous studies concentrated on the standard 
mean field descriptiorP3HH| The interesting point of the Hamiltonian in Eq.|lj) is that under certain conditions on 
the pairing coefficients gk the model is integrable and can be solved exactly with a Bethe ansatz^. This has the 
important consequence that it demonstrates that the BEC-BCS QPT in this model is not an artefact of the mean 
field approximation, even though the ground state correlation functions and energy coincide with the mean field ones 
in the thermodynamic limit for attractive interactions. Furthermore, the defining properties of the Bethe ansatz wave 
function are so qualitatively distinct from its mean field counterpart that it deserves further investigation. Indeed, as 
we will see, some features and insights resulting from the Bethe ansatz state are absent in the mean field solution. 

In this paper we restrict specific calculations to the two-dimensional attractive case where a QPT with unusual 
characteristics signals the transition between two gapped superfiuid phases, one topologically non-trivial, known as 
weak pairing^, and another one characterized by tightly bound quasi-molecules, known as strong pairing^. The 
quantum phase diagram of the two-dimensional p x + ip y spinless fermion pairing Hamiltonian displays several pecu- 
liarities (see Fig.[T]), depending on the fermion density p. For dense gases, p > 1/2, there is no QPT of any sort. When 
p < 1/2, the model displays a QPT between a weak-pairing and a strong-pairing phase, while otherwise the ground 
state energy is an analytic function of the coupling strength g. Volovik anticipated^, using an equivalent effective 
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field theory, that a QPT of this kind can be signaled by a discontinuous change of a topological quantum number. We 
show that the transition is continuous and third order (meaning that the t hird-o rder derivative of the ground state 
energy is discontinuous), contrary to what has been claimed in the literatur e! 21 * 22 *. Furthermore, we demonstrate that 
this phase transition can be related to an experimentally accessible length scale that diverges logarithmically at the 
phase transition point. We give an interpretation of this length scale in terms of the size of the Cooper pairs, with 
strong pairing corresponding to confined pairs that start to deconfine when the system enters the weak-pairing phase, 
and that eventually convert to plane waves in the weakly interacting limit. 




FIG. 1: Quantum phase diagram of the p x +ip y or hyperbolic model in terms of the fermion density p and the attractive coupling 
strength g. It shows two superfluid phases, a strong-pairing (confined) and a weak-pairing (deconfined) phase, separated by 
a third-order confinement-deconfinement QPT at vanishing chemical potential, p, = 0. Cooper pairs deconfine between the 
phase transition line and the Moore-Read line (indicated by a dashed line). The latter line corresponds to p — 1 — 1/ g and 
vanishing total energy, E = 0, and represents a situation where all Cooper pairs are deconfined. The small symbols at p = 0.25 
(quarter-filling) indicate the configurations that are displayed in detail in Figs. [3] and [4] 

A peculiar boundary in the quantum phase diagram is the so-called Moore-Read line^ (shown as a dashed line in 
Fig. [lj. The latter separates in principle a weak-coupling BCS regime from the weak-pairing phase and is defined for 
a particular relation between the density p and the interaction strength g: p = 1 — 1/g. At this line all Cooper pairs 
become deconfined. As explained in detail in the text, our calculations indicate that there is no QPT at the Moore- 
Read line, contrary to what is advocated in Refs. I18|22l Thus, the weak coupling BCS region and the weak-pairing 
phase are adiabatically connected. 

Technically, we solve the more general hyperbolic (or XXZ) Richardson-Gaudin (RG) modeP^ a particular realiza- 
tion of which is the p x + ip y two-dimensional superfluid. We present results for finite systems and, for determining the 
quantum phase diagram, the relevant thermodynamic limit. Calculating those results for the hyperbolic RG model 
implies solving the corresponding Bethe ansatz equations, a task that may happen to be technically challenging, 
specifically for large finite systems. We present a new technique that solves the Bethe ansatz equations numerically in 
polynomial time. For very large systems one can define a thermodynamic limit that allows one to solve the equations 
for the integrable p x + ip y model analytically^]. 

We want to emphasize that any hyperbolic RG model in the continuum limit displays the non-analytic behavior 
observed at the equivalent p, — line as long as the parameters defining the model vanish at a given zero mode (e.g., 
k = mode). The underlying algebraic structure of the hyperbolic RG model allows us to prove a theorem that 
relates subspaccs with different number of pairs but with the same eigenvalues, thereby generalizing a relation from 
Ref. [TS] to all representations of the hyperbolic RG model. These relations are inexistent in the case of the rational 
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RG model (e.g., the s-wave BCS superconductor) and are indeed at the root of the existence of the QPT between the 
weak- and strong-pairing phases. 

In Sec. [IT] we discuss the algebraic properties of the hyperbolic RG model. In Sec. |III| we present numerical and 
analytical results for the p x + ip y model, analyze the phase diagram and interprete the QPT mechanism. Using an 
affinc Lie algebra, we show in Sec. |IV| that the hyperbolic RG model possesses a symmetry that relates the two phases. 
In the Appendices we give a detailed solution of the hyperbolic RG model for one pair (Appendix [A| , for one level 
(Appendix [B]), in the continuum limit (Appendix [C]) and in the mean field approximation (Appendix |D[) . 



II. ALGEBRAIC FORMULATION OF THE INTEGRABLE p x + ip y MODEL AND ITS SOLUTIONS 

A. The generalized Gaudin algebra 

The p x + ip y integrable pairing model can be obtained as a particular parameterization of the exactly solvable 
hyperbolic family derived from the generalized Gaudin algebra. We recapitulate here the main ingredients of the 
procedure to construct such exactly solvable models that we have laid out in previous worlP^. At the same time, this 
serves to demonstrate that particular features of the zero temperature phase diagram of the p x + ip y model are not 
limited to this representation and analogs can be found in many other representations of the same underlying Gaudin 
algebra. Assuming a system characterized by I levels (representations) of an su(2) algebra labeled by the (in principle 
arbitrary) quantum number i, the generic form of the exactly solvable RG Hamiltonians is: 

h = J2 ^1 - E( £ * - e o)x{m, m)sisj -Y.^- e oMm, vj)SfS], (2) 

where the operators S^, S~ and Sf close the su(2) commutator algebra, and X(x,y) and Z(x,y) are antisymmetric 
functions depending on I parameters r\i to be determined later on. The generic Hamiltonian ([2]) commutes with the 
squared spin operator of each level, Sf = S*(Sf — 1) + Sf , and with the total spin operator S = ^ - Let us 
introduce the notation for the spin value of level i, such that (Sf ) = Si(si + 1). Apart from this, one can evaluate 
the commutator between two such Hamiltonians, H and H', based on the same set of parameters r\i but with distinct 
sets of parameters and ej. It turns out that H and H' commute provided that the functions X(x,y) and Z(x,y) 
fulfill the following condition for all x, y and z known as the Gaudin condition^: 

[Z(x, y) - Z(x, z)] X(y, z) - X(x, y)X(x, z) = 0. (3) 

A consequence of this condition is that 

X(x,y) 2 -Z(x,y) 2 =T, (4) 

where T is a constant independent of x and y. Then, by a suitable choice of the parameters in Eq.Q, one can 
define a set of I linearly independent constants of motion, which commute among themselves: 

Ri = Sf - ]T Xfara) {S+SJ + S-S+) - 2 ]T Z^^SfS]. (5) 

3,3^1 3,3¥=i 

Therefore, Hamiltonians of the form of Eq.([2| are integrable^, and furthermore, they can be solved exactly by a 
Bethe ansatz. 

It is convenient to define the following operators: 

S z {x) = -\-Y J Z(x,r lj )S* 1 S±(x) = ^2x(x,r lj )Sf. (6) 

3 3 

They fulfill the commutation relations of an XXZ Gaudin algebra, 

[S z (x),S ± {y)} = ±(X(x 1 y)S ± (x)-Z(x,y)S z (y)), 

[S+(x),S-(y)] - 2X(x 1 y)(S z (x)-S z (y)). (7) 
With this notation, the Bethe ansatz for the cigenstates of order M takes the form 

\*M)=(j[S + {E a yJ\v}. (8) 
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where \v) is a vacuum state such that S i \v) — 0, and S*\v) = —Si\v) for all i. The unknowns E a are solutions of a 
set of nonlinear equations, the RG equations: 

^ Si Z{ m ,E a )- Yl Z{E a ,,E a )= 1 -, Va. (9) 

We will call the variables E a paironP^. Equation Q does not depend on the parameters e^, which means that 
the eigenstates given by Eq.([8| do not depend on these parameters either. This is due to the fact that the ansatz 
from Eq.([8]) represents the eigenstates of the complete set of the integrals of motion and therefore of any linear 
combination of R^s, as given by the Hamiltonian ([2| whose eigenvalues are given by 

E(p M ) = (v\H\v) + 2 SieiZim, E a ). (10) 

i,a 



For r = one recovers the rational RG model, that can be parameterized as 

G 

X(x,y) = Z(x,y) = 

x-y 



(11) 



(any other expression for X and Z that fulfills Eqs.ppt at V = can be transformed to this form through a redefinition 
of parameters). In this work we will focus on the properties of t he hyperbolic RG model, which correspond to T < 0. 
Any set of functions X(x,y) and Z(x,y) that fulfills Eqs.(3p) with T = — 7 2 , can be mapped onto the following 
parameterization^! 



X(x, y) = 27 , Z(x, y) = 7- 



x — y 



x — y 



(12) 



We can simplify the discussion by choosing = Xr/i, with A a parameter that will be fixed later on. Using this 
parameterization, and subtracting a diagonal term 27^ i 77iS|, one obtains an interesting form for the Hamiltonian 
of Kq.lgi: 



H = A(l + 2 7 ) £ mSf ~ A 7 ^(r 7i + m )SfS5 - 2A 7 J2 ^m s t s ] 

i i,3 i,j 

= a (1 + 2 7 - 27^) J2 viS? - 2A 7 sV/A'*, . 



Now we can take advantage of the fact that S z = Y^i <S? is a constant of motion, 

S Z \$ M ) = (M-L/2)\<f M ), 



(13) 



(14) 



where L = 2 Yi Si is the maximum value that the model allows for the order M (therefore L is a useful measure for 
the dimension of the model space; e.g. if Sj = 1/2 for all levels, then L equals the number of levels I). This allows us 
to fix the parameters A and 7 such that 



A (1 + 7(£ - 2M + 2)) = 1, 2A 7 = G, 
leading to a separable Hamiltonian for the hyperbolic RG model: 



The corresponding RG equations, Eqs.(|9|, reduce to: 



1 



7?- — E ■Z—' E , — E 



Q_ 

E n 



= 0, Va, 



(15) 



(16) 



(17) 



with 



(18) 
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Equation ( 10 ) for the energy simplifies to: 

E{<S> M ) = (v\H h \v)+YE a . (19) 



For the remaining discussion we will assume that (y\Hh\v) = 0, which amounts to a simple shift in the energy scale, 
without loss of generality. 

In the pairing representations each su(2) copy is associated with a single particle level i. M is the number of active 
pairs. The vacuum is now defined by a set of seniorities, \v) — \vi, V2, ■ ■ ■ , where the seniority Vi is the number 
of unpaired particles in level i with single particle degeneracy fij, such that Si — (51^ — 2^j)/4. 

In two spatial dimensions, one can define a representation of the su(2) algebra in terms of spinless fermions in 
momentum space, c k , c k . One obtains a level r\k for each pair of states (k, —k), where the index k now refers to the 
momentum in two dimensions (in order to avoid double counting we select k x > to label the levels). Furthermore, 
one can include a phase factor in the definition of S^: 

q Z l/'t it l\ Q+ ^ x V t t C— x V IO(W 

2 \ feCfe C -k c -k — i J > h k — Hu c fc c -fe' ^fe — Hu c -k c k- ( Z[) ) 

By taking rjk = k 2 , one obtains the p x + ip y model presented by Ibanez et. aP^. 

A- 2 



H Pa .+i Py 



y ( c l c k + c -k c -k) ~G ( fc - + ik v)^ ~ lk 'y) c l c -k c -k' c k>- (21) 



fe,fcx>0 fc,fe :c >0, 

fe',fe' >0 



In the case of cold atom gases, the spinless fermions operators represent atoms in the same hyperfine state. 

B. Singularities in the Richardson-Gaudin equations 

J ust lik e in the rational case, the RG equations become singular when two or more pairons approach the same 
leve l 27 ! 2 ^ . Analyzing the residues in the equation, one finds that the number of singular pairons has to be equal to 
2sfc + 1 for a singularity around level 7%. However, around E a — the Eq.pT| for the hyperbolic RG model displays a 
peculiar singularity not present in the rational models. Let us analyze what happens when N of the pairons converge 
to zero energy, and M — N pairons do not. For the pairons that do not converge to zero, one finds that the equations 
correspond to a state of order M — N with all pairons distinct from zero. In the equations for the N singular pairons 
the contributions of the non-singular pairons and the levels are of order one and become negligible compared to the 
singular terms. Therefore we can replace the non-singular terms by a single level, with an averaged level parameter 
i] given by 

Appendix [B] gives the analytical solution of the one-level problem. One finds that a necessary condition to have N 
values converging to zero is that 2Q + 1 = N, or 



1 

G 



= L-2M + N + 1. (23) 



We recognize two special cases: 



special case (i) N = 0: 

no pairons can converge to zero for q < L — 2M + 1; we will see that this boundary coincides with the phase 
transition line in Fig. [T] 

special case (ii) N = M; 

all pairons can converge to zero for g = L — M + 1; this situation determines the so called Moore-Read lin e! 18 * 23 ! 
in Fig. |T] 

Between these two regimes, a fraction of the pairons can converge to zero at integer values of G _1 . 

This brings us to a qualitative view of the quantum phase diagram of Fig. [TJ on the left side of the diagram one 
finds the weak coupling BCS region, where the pairons behave much like in the rational model: they form arcs in the 
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complex plane. If the coupling strength is increased, one reaches the Moore- Read line, that coincides with special case 
(ii) where all pairons collapse to zero. Increasing the interaction strength, still in weak pairing, some of the pairons 
can collapse to zero (at integer values of G _1 ). Finally, for strong interactions and M < L/2, one can pass special 
case (i) and make the transition to the strong-pairing phase. Here the ground state corresponds to a configuration 
where all pairons are real and negative. 

At the Moore-Read line all pairons E a converge to zero, which means that the total energy E also goes to zero. 
This coincides exactly with the mean field result. The corresponding wave functions are not the same because the 
mean field wave function breaks pair number symmetry, and can be written as a superposition of exact zero-energy 
states for a range of numbers of pairs: 



|V>mf) oc exp J A T~~]r c i c - 
\ k k >o x v 



k (24) 



while the exact wave function is given by the projection of \ip m f) onto the M-pair sector, 

M 



\k,kx>0 " 



c{cl k \ W- (25) 



Just below the line, when N = M — 1 = L — 1/G — 1, one finds that the energy of the exact eigenstate is equal to 
the energy E p of the one-pair state calculated in Appendix A~| because of the symmetry to be shown in section IV 



Increasing the number of pairs by one while maintaining G fixed, moves the state to the Moore-Read line where the 
total energy is zero. This means that for a small increase in density, AM/L = l/L (which becomes infinitesimal in 
the thermodynamic limit), there is a finite jump in energy, AE = \E p \. This is interpreted by Ibahez et al. as a 
signature of a zero th order QPT. Note however, that E scales proportional to wL, where w is the kinetic energy cutoff 
as defined in App. |Aj Therefore, the shift in energy AE resulting from a shift in density AM/L is given by 

which is of order uj for AM = 1 because the derivative is free of scales. This means that a jump in energy of size \E p \ 
for a change in density of size 1 / L is not a sign of a phase transition but rather a consequence of the pathological 
scaling behavior of the energy in this model. 

III. QUANTUM PHASE DIAGRAM OF THE INTEGRABLE Px + i Py MODEL 

A. The case of finite systems 

We would like to stress here that in order to obtain the exact solution of the pairing eigenproblem for M pairs in a 
system of size L, one has to diagonalize a Hamiltonian matrix of dimension L\/(M\(M — L)\). The integrability of the 
RG model reduces the exponential complexity of the problem to that of solving a set of M non-linear equations in M 



unknowns, Eq.( 17), which can be done in polynomial time. This justifies labeling these models as exactly solvable. For 
instance, below we show exact results for L = 504 and M = 126, which corresponds to solving an eigenvalue problem 
of dimension 10 122 , while in our case each eigenvector arises as a particular solution of the 126 coupled non-linear RG 
equations. 



Although easily written down, actually solving Eqs.(17l for more than a few pairs turns out to be a cumbersome 



task due to singularities that occur when some of the pairons approach the origin or a certain level value rjk, or if 
two or more pairons approach each other. For the rational model with all = 1/2, wh ere the singularities entail at 



most two pairons, Richardson already proposed a practical solution 2 ^, by rewriting Eqs.(17) in terms of the real part 
of the pairon values and the square of the difference between two pairon values. If singularities with three or more 
pairons occur, this trick no longer works . More complicated changes of variables have been proposed to remediate 
this problem for the rational modePZES, but in the hyperbolic model we are confronted with a situation where all 
pairons can converge to zero and where the change of variables can no longer be performed with sufficient accuracy. 



Standard gradient methods to solve the Eqs.(17) directly run into problems because the Jacobian matrix becomes 
ill-conditioned well before the pairons reach the singularities. 

Our numerical strategy to solve the RG equations is to start from an educated guess for the position of the M 
pairons in the weak coupling limit. In this limit the RG equations decouple into M independent equations (one for 
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0.04 


0.08 


0.16 


0.20 


0.32 


0.36 


0.40 


0.52 


0.64 


0.68 


0.72 


0.80 


1.00 


Sk 


1 


1 


1 


2 


1 


1 


2 


2 


1 


2 


1 


2 


3 



TABLE I: Level parameters r/k and Sk for a disk with a radius of five unit cells in a two-dimensional square lattice. 



each level). As discussed in Appendix [Bj the solution for each of these generalized Stieltjes equations is given by the 
roots of a Jacobi polynomial, which we use as a starting point for the iterative procedure at small values of G. The 



value of the coupling strength G is then gradually increased up to the desired value. At each step we solve the Eqs.( 17) 
using a variation of the Levenberg-Marquardt algorithm^], where care has been taken to avoid singularities in the 
Jacobian matrix. In order to avoid the typical divergencies that burden the numerical procedures in intermediate 
steps, we modify the level parameters rj k and the central charge Q by adding a modulated artificial imaginary part. 
Once the desired value of G is obtained, a new iterative procedure is initiated, in which the imaginary parts of rj k 
and Q are gradually reduced to obtain the final solution. For the numerical results displayed in Figs. [2]to[TlJ we have 
taken the kinetic energy cutoff ui = 1/2, such that the level parameters rjk range from to 1. 




FIG. 2: Real and imaginary parts of the spectral parameters E a as a function of g = GL for 10 pairons in a disk of dimension 
L — 40 (quarter filling), with levels as specified in Table [I] 



As a first illustration of the properties of the exact solution in the different regions of the quantum phase diagram, 
we treat ten pairs in a two dimensional lattice in a disk geometry with a five unit cell radius. The system has a 
total active space of dimension L = 2^ fc s k = 40, hence this amounts to quarter filling, i.e. p = 0.25. The resulting 
level parameters are depicted in Table |TJ In the extreme weak coupling limit the ten pairs fill the lowest four levels. 
With increasing coupling the pairons expand in the complex plane lowering their real part, and consequently the 
total energy. There is a particular value of the scaled coupling strength g = GL at which all pairons collapse to zero 
(the Moore- Read line), later on expelling one pairon with a negative real value. From this point up to the QPT at 
Q = there are a series of collapses with nine, eight, etc. pairons, and after each one a new real negative pairon 
is produced. The exact solution for this small system shows the qualitative behavior of the pairons in the quantum 
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phase diagram (Fig. [I]). In order to get a more quantitative picture of the pairon distribution in the three regions of 
the quantum phase diagram, we plot in Fig. [3] the pairon distributions for three representative values of the coupling 
strength, g = 0.5, 1.5, 2.5, at quarter filling for a disk of radius 18 corresponding to a total pair degeneracy L = 504, 
indicated in Fig. [I] by their respective symbols- 1 . One notes the distinctive features: in weak coupling BCS part of the 
pairons stick to the lower part of the real positive axis, while the remaining pairons form a double arc. Approaching 
the Moore-Read line it looks like the arc is going to close around the origin, but just at the Moore-Read line all the 
pairons collapse to zero. A first real negative pairon emerges. In the intermediate weak-pairing region a successive 
series of collapses ensues, at integer values of Q, each time producing one more real negative pairon and reducing the 
size of the arcs around the origin. When the last pairon turns real and negative, the system enters the strong-pairing 
phase. From then on the most negative pairon diverges proportional to the interaction strength G, while the least 
negative pairon converges to a finite value that can be related to the condensate fraction, as we will see later on. 
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FIG. 3: Pairon distribution for L = 504 at quarter filling, p — 0.25, for g = 0.5, g = 1.5, and g = 2.5. 

In Fig. [4] we show the momentum density distributions corresponding to each of the configurations in Fig. [3j In 
the weak-pairing region one sees that the lowest momentum states are fully occupied. On the other side of the 
phase transition the lowest momentum states are vacated, resulting in a non-trivial topology for the two-dimensional 
momentum distribution. 
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FIG. 4: Momentum density profiles, n*,, for the same parameter sets as in Fig. [3] (see symbols), compared to the mean field 
solution for the same system (continuous lines). 



B. The thermodynamic limit 

We are interested in determining the quantum phase diagram of the integrable p x + ip y model defined by the 
Hamiltonian ( |21[ ). To explore the possible QPTs we have to properly define the thermodynamic limit for this model. 
From Fig. [5] we see that for large systems the exact solution and the mean field solution tend to the same results. 
Particularly, in mean field theory one observes that for static values of the chemical potential /i and the pairing field 
A, the quantities M, E and G _1 scale proportional with L. We can describe the thermodynamic limit when L goes to 
infinity in terms of intensive quantities such as the scaled interaction strength g — GL and the density p — M/L. At 
the same time we have to establish a cutoff u> in kinetic energy, such that the energy density e — E/L remains finite, 
and a density of levels g(rj), such that g{rf)drj = 1. The cutoff u> determines the energy scale of the problem. It 
can be renormalized by relating it to the energy of a single pair as explained in Appendix [A] 

As demonstrated in Appendices [C] and |Dj the exact solution and the mean field solution coincide in the thermody- 
namic limit and yield the following integral equations: 

*%(»/) \ ===dn = -, (27) 

2w Ou 1 

Q(V) „ „ = =m d ^ = ~-l + 2p. (28) 



o 



o 



Note that the chemical potential fi and the pairing field A 2 scale proportional to uj, and can be determined from the 
self-consistent solution of the above equations (see Figs [6] and [7| . The energy density in the thermodynamic limit is 
given by 



uj A 2 1 /* 2lj i 

e = - + M2p - 1) + / g(v)\/(v - 2^) 2 + 4 ?7 A 2 d7 ? . (29) 



9 2 

One sees that Eqs.( 27|28 ) coincide with the conditions 



(i 



de de , 



meaning that the mean field solution minimizes the energy written as a function of the mean field parameters /z and 
A. 
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FIG. 5: Finite-size scaling of the energy density e at quarter filling for several interaction strengths (corresponding to distinct 
regions of the phase diagram), both for the exact solution (symbols) as for the mean field solution (lines). 




For the two-dimensional p x + ip y model one finds that g(rj) = l/(2w), and one can solve the integral appearing in 
Eq.(29| analytically: 

w A 2 |/i|(A 2 -/i) w + A 2 -u n -7, t^t 

+ A2 ( A2 -^lnf- + A2 -^^V;) 2 + 2 - A -V (31) 
2oj I A 2 -^+|/x| J x ' 
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FIG. 7: Analytical results for the pairing field A and its derivative dA/dg in the inset, as a function of g for various densities. 



From Eqs.(30) we obtain explicit expressions for g and p: 
1 



.9 



P = 



1 1 

2 _ 2lu 



v/(w - p) 2 + 2uA 2 - \n\ + (p - A 2 ) In 



- p) 2 + 2uA 2 - \p\ - A 2 In 



- p + A 2 + y/(uj - p) 2 + 2cjA 2 

- ji + A 2 + y/(cj - p) 2 + 2^^ 

a 2 -m + H ; 



This allows us to rewrite Eq.(31 1 in a more elegant form 



where 



. A 2 ujA 2 m+I^I 
2g fi 2 



(32) 
(33) 

(34) 
(35) 



The appearance of the absolute value of \i clearly indicates that there exists a non-analyticity in e at [i = 0. 



C. Derivatives of the energy density 



The third term of Eq. ( 34 1 might suggest a singularity in e at fj, 

coincides with q 



and hence a zeroth order QPT at this point. 
such that the singularity in e is canceled out and as a result e 



However, it turns out that /i 

is a continuous function of g and p over the whole parameter range. A special point that one observes in Eqs.( 32p3 ) 



corresponds to p = A /2, which leads to g = 1 — p. This is the Moore-Read line we discussed in Sec. II B As argued 
above and as can be seen from Fig. [HJ the energy vanishes at this point but i s cont inuous. Consequently, there is no 
zeroth order QPT in this model (contrary to previous reports in the literatur e! 18 * 22 ^). 

The fourth term in Eq.(34| might suggest a discontinuity in the derivative de/dp at p = 0. However, it turns out 
that the third term also has a discontinuity in its derivative, that exactly cancels the one in the fourth term. To 
formalize this, let us consider the derivatives of e with respect to p and g: 



de 
dp 



2/i, 



de 
dg 



= -A 2 / 5 2 . 



(36) 



In Eqs.( |32|33 l we see that g 1 and p are continuous functions of p and A. Hence, p and A are continuous functions 
of p and g, and therefore the derivatives of e with respect to p and g are continuous functions of p and g. We conclude 
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that there is no first order QPT in this model. Note also that in the strong coupling limit A tends to the same value 
for p = 0.25 and p = 0.75 (see Fig. [7]). This is a consequence of the symmetry that we will discuss in Sec. IV 

For the higher order derivatives, we resort to a graphical representation. From Fig. [9] we see that th ere is a 
discontinuity in the third-order derivative of e. The discontinuity occurs when p = 0. If we set p = in Eqs.( 32p3 1, 
we find that th is coi ncides with = 1 — 2p. This transition line has been identified in previous works as the 
Read-Green lin e! 14 * 18 !. We see that a third-order QPT occurs at this line, for densities below half filling, i.e. p < 1/2. 




FIG. 9: Higher-order derivatives of the energy density e as a function of g for various densities. The small circles mark the 
transition point at g = (1 — 2p) _1 . 
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D. Characterization of the superfiuid phases 



In s-wave pairing models (rational RG model) the transition from a weak to a strong-pairing phase is often described 
in terms of a "BCS-BEC" crossover, where on the weak side one has dilute pairs moving in a coordinated way as 
described by BCS theory, while on the strong side the fermion pairs behave as molecules that become perfect bosons 
in the strong coupling limit, where they form a pure BEC. In the p x + ip y case the picture differs in an essential way: 
we see that instead of a crossover there is a sharp third-order QPT at g^ 1 = 1 — 2p. 

Because the exact ground state correlation functions coincide with the mean field ones in the thermodynamic limit, 



we can use Eq.(D17) to evaluate the condensate fraction n c , following the standard definition of Yan 



where a and b are such that 



fe,fc x >0 

LA 2 f 2 " r, 
2uMj (r, - 2/i) 2 + 477A 2 V 
A 2 oln(l - 2w/a) - 61n(l - 2u/b) 
2up a — b 



(?y - 2pf + ArjA 2 = (r]-a)(r]-b). 



(37) 



As explained in Appendix [C] a and b are the lower and upper bounds of the line of pairons in strong pairing in the 
thermodynamic limit. We have used the mean field expectation values because they coincide with the exact results 
in the thermodynamic limit. For finite systems one could evaluate the exact correlation functions using techniques 
similar to the ones presented in Ref. |33]for the rational RG model. Figure [TU] displays the results for various densities. 
One observes that a pair condensate develops at large interaction strength. Two features draw attention: first of 
all, we see that for a fixed density p < 1/2 the condensate fraction reaches a maximum in the strong-pairing phase 
at a finite value of g, and then slightly decreases to its limiting value for g — > 00. This means that the fermion 
pairs do not behave as pure bosons, but rather as bound molecules with a certain spatial extension, such that the 
Pauli principle prohibits a complete condensation of all molecules^. Secondly, note that because we study here the 
condensate fraction as a function of the interaction strength, the limiting value of n c in the strong BEC limit is not 1, 
contrary to what one is accustomed to in three dimensions for the s-wave case^, where one evaluates n c as a function 
of the ratio between interparticle distance and scattering length. 
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FIG. 10: Condensate fraction of pairs, n c , as a function of g, for densities p = 0.05, 0.10, . . . , 0.95, in descending order. The 
circles mark the transition points at g — (1 — 2p)~ 1 , while the triangular symbols mark the point where n c is maximal. 

To understand the behavior of the condensate fraction in the strong coupling limit, we evaluate p. and A 2 for j> 1 
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up to leading order in g: 



p 2ujg{\ - 2p/3)(l - 2p)/4, (38) 

A 2 - 2ug 2 {l - 2p/3) 2 p/2. (39) 

Using the expressions for a and b in terms of p and A 2 from Appendix [cj and substituting these in Eq.(37l, one finds 
that in the strong coupling limit 



lim n c 

L— >oo 1 g— ¥ oo 



1 

4p 



l-il^!ln 
Ap 



l + 2p 
1 - 2p 



(40) 



At low densities, this expression approaches Yang's bound for the condensate fraction of a generic fermionic 
superfluicP^, n c < 1 — p. 

Given that there is a condensate, it is worthwhile to analyze the corresponding condensate wave function in mo- 
mentum space, which in the thermodynamic limit is given by^l 



^Yang(fc) = U k V k OC 



h x thy 



^/(fc 2 -a)(k 2 -b)' 

The condensate wave function has an inherent length scale r rms that can be calculated ad^H 

\V(f){k)\ 2 dk/ J \cj)(k)\ 2 dk. 
After some algebra, one obtains the analytical expression: 



2lo [ab — (a + 6)cj] 



(2ui-a) 2 



(2cj-b) 2 



a ~~a- b b +b I ln C 1 - 2tjJ /a) - In (1 - 2u/b)] 



rms. Yang 



(a - b) [aln(l - 2u)/a) - 61n(l - 2w/b)] 



(41) 



(42) 



(43) 



The resulting radius r rmSj Yang is plotted in Fig. 11 One sees that it goes to a finite value in the strong coupling limit, 
and that it diverges at the QPT, with a singularity in the function ln(l — 2u>/a), that scales as In \p\. The condensate 
wave function can be related directly to observable quantities: from Wick's theorem for the mean field solution one 
obtains 



4 C r'/ ( C r' C r) = (4 C r' C r' C r) ~ ( c £c„ 



C l' C r') + (4<V 



We can rewrite this as 



^Yang^ - r')\ 2 oc (^(n r - (n r )) (n r > - (n r /))^ + \F(r - r')\ 2 , 



(44) 



(45) 



where 0Yang(^) is the condensate wave function in coordinate space and F (r) is the Fourier transform of the momentum 
density \c\c k \. In trapped Fermi gases the first term of the right hand side of Eq.(45) can be obtained through 



quantum noise interferometry251, while the momentum density and hence F(r) can be obtained from a time-of-flight 
analysis after releasing the trap. Once these quantities are known, one can calculate the radius r" rmSi Yang- Therefore 
this length scale is an experimentally accessible quantity that would give a clear signature of the third-order QPT in 
the p x + ip y superfluid. 

We can apply a similar analysis to the Bethe ansatz solution: in Eq. ^ we see that the exact eigenstates are given 
by a product of pair wave functions, 



S+(E a ) = ^0 BQ (fc)c t fe ci fe , with Mk) 



k 2 -E 



(46) 



The ground s tate is given by a set of pairon values E a and hence can be interpreted as an ensemble of Cooper pairs 
of varying siz e 1 35 1 37 1 . Evaluating the root mean square radius for the Cooper pair wave function </>_e(r) in coordinate 
space, one finds that 



2wRe [E] (Re [E] - 2w) f ^ 



\E\ 2 [ln(l - 2lo/E*) - In (1 - 2lo/E)} 



rms, E 



(E - E*) 2 [E\n (1 - 2lu/E) - E* In (1 - 2u/E*)] 
4uj (3E 2 - QujE + 8uj 2 ) 



3E(E ~ 2uf [ln(l - 2u/E) + 2u/(E - 2u)] ' 



for E real and negative. 



(47) 
(48) 
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This length scale is finite if E is complex, and converges to a finite value when E approaches the real negative axis. 
But when approaching real and positive values of E, inside the interval [0, 2w], the length scale r rm5j e diverges. The 
exact solution of the integrable p x +ip y model gives us a beautiful insight into the nature of the QPT: in strong pairing 
all pairon values E a are real and negative, which means that the exact ground state corresponds to an ensemble of 
confined Cooper pairs with finite radius. Lowering the interaction strength one encounters the QPT at which one 
pairon value becomes real and positive, which corresponds to a single deconfined Cooper pair on top of the ensemble 
of confined Cooper pairs. Entering further into the weak-pairing phase, the pairon values become complex. At 
integer values of Q there are 2Q + 1 deconfined Cooper pairs, while at intermediate values part of the pairons form 
two arcs around the origin in the complex plane, some stick to the lowest part of the real positive axis (inside the 
arcs) and the remaining ones stay confined on the real negative axis. At the Moore-Read line all Cooper pairs are 
deconfined. From there on part of the Cooper pairs will form resonances of finite size (complex pairon values), till all 
Cooper pairs convert to plane waves in the extreme weakly interacting limit. Hence we can understand the QPT as 
a confinement-deconfincmcnt phase transition. 



10 i- 
8 - 




10 



9 



FIG. 11: Root mean square radius, r rms , 
lines indicate the Moore-Read point (e 



of the condensate wave function in units of 1/ 
0), and the phase transition point (/x = 0). 



2w, at quarter filling. The dashed 



IV. SYMMETRIES OF THE RICHARDSON-GAUDIN EQUATIONS 

In this section we highlight some symmetry relations of the RG equations that allow us to relate states in the 
strong-pairing phase with mirror states in the weak-pairing range that have the same energy^. To set the stage we 
would like to introduce a new set of operators, 

{T+,T-,T*}, with m = ■■■ ,-2,-1, 0,1,2,.-. , (49) 

which form an affine Lie algebra 

A possible representation of these elements of the algebra in terms of su(2) spin operators is (a = z, ±) 

is = £fa) m/2 sf , (si) 

i 

for (in principle arbitrary) quantum number i. 

It is straightforward to check that the following two families of Hamiltonians 

H r = Ti m -GT+Tu, (52) 
Hh = T2 m — GT+T m , (53) 
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define the integrable rational and hyperbolic RG models of Eq.( 16 1. It is illuminating to solve those integrable rational 
and hyperbolic models in terms of the generators of the affine Lie algebra. The eigenvectors can generically be written 
as in Eq. |8| with (unnormalized) 



S ± (E a ) 
S ± (E a ) 



G 



T T 2m i E a ^ , (0) = T% , rational model, 



E 

m=0 

E T 2™+i ' E « * ' S± (°) = T -i ' hyperbolic model, 



m=0 



17 11 < 



(54) 
(55) 



where E a are the pairons. Notice that to represent these eigenvectors one needs an infinite set of generators, the 
even series for the rational model and the odd for the hyperbolic model. In terms of the affine operators one can 
demonstrate some interesting symmetry properties of the model: 

Theorem 1: Given the hyperbolic Hamiltonian Hh, and a particular eigenstate |$m) of order M < (L — G~ 1 )/2, 

for integer values of G _1 = L — 2M + 1 — A, there exists another eigenstate \$>m+n) = |$m) with the same 

energy E. 

Proof: Define Hh = T 2 Z — G T^T^ . One can straightforwardly evaluate the commutators 



[i^T+j = T++2GT+T$, 
[[ff fc ,T+] ,T+,] = 2GT+T+,. 



Now let us consider a state of the form 



1$ 



M+N/ 



MI 



(56) 
(57) 

(58) 



We want to determine the conditions for \&m+n) to be an eigenstate if |$m) is another eigenstate of Hh of energy 
E. Using the commutation relations above we have 



H h ,{T+ 1 ) N =N(T+ 1 )" 'T+[1 + 2GT Z +G(N-1)} 



, N-l 



Therefore 

H h \$m+n) =E\§ m +n) 
The condition to be an eigenstate is 



, JV-l 



=L-N-2M+1. 



(59) 



(60) 



(61) 



This means that for a certain value of G = 1/(L — M — M' + 1), there exist two states, one of order M and another 
state of order M', that have the same energy E. From this algebraic relation, one can distinguish the two special 



cases already mentioned in Sec. II B 



Special case (i) A first case is the symmetric case, that occurs when M — M', i.e. when 

Af = (L + l-G- 1 )/2, Q = -i; 



(62) 



Special case (ii) Another special case occurs when M' = 0, which amounts to 



M = L + 1- G"\ 



M — 1 



(63) 



, M 



The corresponding eigenstate of order M is given by \0m) = (T^Y" \v). The energy of this state is equal to 
the vacuum energy, E(0m) = (v\Hh\v). Comparing this state with the Bcthc ansatz from one observes 

that this state corresponds to the case where the M pairons converge to 0. 



These algebraic relations are easily recognized to be symmetries of the RG equations as well. Consider Eqs. (17) 
and Q for G^ 1 = L - 2M + 1 - N, and evaluate Q and {l8} for the cases with M and M + N pairs, where N 
pairons Ep, = M + 1, • • • , M + N, are zero. Then, Q for the M pair case is given by Q = —(A + l)/2 while for 
the M + N situation it is Q = -(N + l)/2 + N. On the other hand, the first term in Eq. (TP?) is the same for M 
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or M + N but the second one differs precisely by N/E a , so that it cancels with the —N/E a coming from the third 
term. The bottom line is that the RG equations remain the same for the two cases M and M + N. 

In the thermodynamic limit this symmetry mani fest s itsel f as a particular transformation of the chemical potential 
\x and gap A which preserves the form of the Eqs.(27) and (28). For certain values of /j, and A 2 there exists another 
set of values // and A' 2 that lead to the same equations, 



(j, 
A' 2 



A 2 - 



2(i. 



(64) 
(65) 



The theorem above relates eigenstates belonging to different Hilbert subspaces for particular values of G and is a 
result also obtained in Ref . in the context of the p x + ip y superfluid. By using the affine Lie algebra introduced in 
Eqs. (50), we have generalized this result to arbitrary realizations and representations of the general hyperbolic RG 
model. One may wonder whether a similar symmetry property holds for the rational RG model H r = 



One may wonder whether a similar symmetry property holds for the rational RG model H r 
Let us assume that there exists a state 



GT+T -. 



'M+N) = {-L-2) I^M/ , 

and determine the conditions to be an eigenstate of H r whenever |$m) is- The commutator 



[H r ,(T^) N ]=N(T+ 2 f- 1 [H r ,T± 2 ] 



N(N 



il(T+) Ar - 2 [[S r ,T+],T+ 



(66) 



(67) 



implies that the action of the Hamiltonian on the state ( 66 ) is 



H r \$m+n) = E \$> m+ n) + N (T+ 2 ) N 2 T+ [T+ 2 (l + 2GT* 2 ) + G (N - 1) T+ ] |$ 



Ml 



(68) 



To cancel out the second term two conditions need to be satisfied. First, the number of pairons with E a = must 
be TV = 1. Second, 



'l + 2GT1 2 ) |$ M > = 0, 



(69) 



which indicates that for a particular value of G, satisfying the relation above, a state with at most one zero pairon may 
exist in the rational model, as opposed to the hyperbolic case where states with an arbitrary number of zero pairons 
are possible eigenvectors. This qualitative difference reflects in the quantum phase diagram. While the hyperbolic 
model displays a QPT from a weak-pairing phase to a strong-pairing phase, the rational model displays a crossover. 
This is so, since the symmetry relation indicates that at the /i = line the ground state becomes degenerate, while 
this is not the case in the rational model. Connecting points with equal energies in the phase diagram, we see that 
they form two lines (one in weak pairing and one in strong pairing) that approach the phase transition line at the 
same point but under different angles. Hence it is no surprise that the energy is continuous but has a discontinuity 
in the third and higher derivatives at this point. 

Another interesting symmetry relation, not discussed before, emerges by considering the integrable Hamiltonian 



H h = Tf 



G T+i T 



15 



(70) 

A simple way 



which is another exactly solvable hyperbolic RG Hamiltonian, corresponding to parameters fji = 1/rji 
to derive it is to make the choice e, = A?^ 1 in Eq.([2]). 

Theorem 2: Given the hyperbolic Hamiltonian Hh, there exists another Hh with G _1 = L — 2A1 + 2 — G _1 whose 
eigenstates are the same but the eigenvalues are given by E — (v\Hh\v) + J2 a 1 - Hence, although Hh must have 
the same eigenstates as Hh for the same Hilbert subspace of order M, their eigenvalues are different. 

Proof: Wc explore now how the eigenvalues of Hh are related to those of Hh- The commutator with Hh is given by 



H h , H h \ = {T+TZ 1 - T+ Tf) [2GG (1 — T Z ) —G — G 
Acting on a state of order M it results 



Hh, Hh 



\<l> 



Ml 



GG (L — 2M + 2) - G — G (T+TZ 1 - T+ x Tf ) |$ M ) 



Therefore, for given values of G and M the relation 

G" 1 + G" 1 = L - 2M + 2, 



(71) 



(72) 



(73) 
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defines a value of G such that [Hh,Hh] = 0. This algebraic identity also shows up in Eq.(17), since these equations 
are invariant under the transformation 

fji = %\ E a = E-\ G- 1 =L-2M + 2-G-\ (74) 

Hence one can conclude that for a given eigenstate \$m) of Hh, the following relation holds: 

a 

In particular, at G _1 = L — 2M + 2, the state |$m) becomes an eigenstate of the operator T^ 1 TZ 1 - 

V. CONCLUSIONS 

The exactly-solvable hyperbolic RG model allows a parameterization that leads to a pairing model for spinless 
fermions in two dimensions with a p-wave interaction of p x + ip y character. The model can be solved with a Bethe 
ansatz, even for large system sizes, using a new numerical technique presented here. The thermodynamic limit can 
be evaluated analytically, leading to expressions for the ground state energy and few-body correlation functions that 
coincide with mean field theory. 

Even though the p x + ip y model has been studied before^^H, we reach drastically different conclusions about 
the zero temperature quantum phase diagram. Our numerical and analytical results show that the quantum phase 
diagram exhibits a third-order quantum phase transition when the chemical potential changes sign. Most importantly, 
as shown in this paper, this transition is characterized by a logarithmically diverging length scale that can be obtained 
experimentally from correlations in the density fluctuations. This length scale can be associated with the size of the 
Cooper pairs making up the correlated many-body state. Those Cooper pairs, tightly bound in the strong-pairing 
phase, start to deconfine at the QPT point where /i = 0. Thus, in a certain sense, our physical picture is that the 
strong-pairing region represents a confined superfluid phase while the weak-pairing region is the deconfined phase. 
Similarly to many other examples of phase transitions which are not of the Landau type, as for instance in lattice 
gauge theories, the transition is not associated with any broken symmetry and thus has no local order parameter. 

We show that these are general features of the hyperbolic RG model, irrespective of the specific symmetry or 
dimensionality of the representation. This means that our conclusions equally apply to other models derived from 
the same generalized Gaudin algebra, be it models for pairing in atomic nuclei, or Jaynes-Cummings like models of 
resonantly coupled fermionic atoms. Furthermore, the algebraic underpinnings of the hyperbolic RG model allow us 
to demonstrate that there exists a transformation that relates states in strong pairing with states in weak pairing, 
and that the phase transition line corresponds to the states that are symmetric under this transformation. 

Finally, we would like to emphasize the fact that the phase diagram is based on the exact solution, which guarantees 
that the phase transition is not an artefact of the mean field approximation. Moreover, the Hamiltonian ([!]) is 
integrable for attractive and repulsive interactions while there is no mean field solution in the yet unexplored repulsive 
case. Therefore, exact solvability could be a unique tool to investigate repulsive p-wave interactions. 
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Appendix A: One pair 



Already the state with one pair leads to a non-trivial result. Equation (17) reduces to: 



Ei 2 ^ Ei - m. 

In order to solve this equation, we have to specify the space in which the model is defined. Let us consider a disk in 
momentum space, with radius aR such that R is integer and a is the spacing between successive momentum states. 
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The kinetic energy dispersion is given by k 2 /2, which means that there is a kinetic energy cutoff lu — (aR) 2 /2. In 
the continuum limit, corresponding to R 3> 1, one can take L — irR 2 /2 and replace the sums over k by integrals with 

clksg — dky — Qj . 



E /( fc2 



f(k 2 )dk x dk y 



k x >0,k<aR 



all 



f{k 2 )kdk 



k=0 



irR 2 /2 [W) 



(aR) 2 Jo 



f(r))dr) 



2cu 



2uj 



(A2) 



Hence, 



E 

k,k^>0 



1 



L 



Ei - k 2 2w J Ei-r) 



1 



-dr\ 



L f 2oj\ 

In 1 

2uj \ Ei J 



Substituting this result in Eq.(Al|, one obtains that 



^'I-K 1 -!) 



o. 



In the limit L> 1 one finds that E\ = 2uj/ f g , where g = GL and f g is the solution of the following equation: 

(l/ 9 -l)/ 9 = m(l-/ g ). 
If l/g < 1, then f g <0 and one obtains a ground state energy for one pair, 

E 1 = 2u/f g , 



(A3) 

(A4) 

(A5) 
(A6) 



such that Ei < 0. This value Ei then defines a physical energy scale that can be used to eliminate the cutoff ui as 
done in Ref. [THl For l/g > 1, the one-pair ground state corresponds to a non-paired state with Ei = 0. In this case 
there is no relevant energy scale to renormalize the cutoff. 

The one-pair energy is a physically meaningful quantity, that allows one to relate different models describing the 
same physics. Typically, models with p x + ip y pairing will have an interaction that scales as (k x + ik y ){k' x — ik' y ) 
for small fc, just like the integrable p x + ip y model, but for large momenta the interaction could be very different. 
Because the peculiarities of the quantum phase diagram are determined mainly by the small-fc behavior, we expect 
similar results for other two-dimensional p x + ip y models, as e.g. in Ref. 1161 To relate their results to ours, one has 
to renormalize the interaction strength, which can be done by equating the one pair energy. 



Appendix B: Polynomial solution for the one-level model 



In the limiting case of a single level, I = 1, Eq.(17) takes the form: 

_J v 1 = Q Va 



(Bl) 



Equation (Bl I is a generalized Stieltjes equation^, that can be solved by multiplying it by 2E a (rj — E a ) P' M (E a ), 



where Pm{ x ) — Ila^i^ ~ E a >) is the polynomial that has the unknowns E a as its roots. One finds that 

[2(a + Q)E a - 2Qt 1 ] P' M {E a ) + E a {r, - E a ) P^{E a ) = 0, (B2) 
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for a = 1, . . . , M. Since this is a polynomial expression of order M, valid in M distinct points (assuming that none 
of the E a coincide), one can identify the polynomial expression with a multiple of P{x): 

[2(s + Q)x - 2Q V ] P' M {x) +x(j]-x) P^(x) = M [2(s + Q)-M + 1] P M (x), Vs. (B3) 

This equation can be transformed into the Jacobi equation in the variable y where x = (t?/2)(1 + y). The solution is 
a Jacobi polynomial: 

P M (x)ocV^(2x/r)-l), (B4) 

with parameters a = 2Q/r; — 2Q — 2s — 1, /? = —2Qjr\ — 1. 

Let us expand Pm(x) = ^2iLoPM,iX l , where pm,m = 1- By equating orders of a; z in Eq.(B3), one obtains a 
recurrence relation that defines the coefficients Pu,i completely: 

PM,i-i = V [M + i _ 2 ~Q)]{M + 1 - if M ^ ^i = M,M (B5) 

If 2Q + 1 equals an integer value N < M, one sees that Pm{x) has a factor x N . This means that a singularity occurs, 
in which N pairons will collapse at 0. Of particular interest are the cases N = M, were all pairons collapse to zero, 
and N = 0, where no more singularities occur. 

For multilevel systems, e.g. a system of I levels, the polynomial Pm(x) is still well defined. Rather than a Jacobi- like 
equation, it now has to fulfill a generalized Lame differential equation with negative residues^, 

P 'm{*) + (E -\] PM ~ n / A ; (X) Pm{x) = 0, (B6) 

with ro = —2Q, r^ = —2sk, and 770 = 0. The solutions Pm(x) are known as Stieltjes polynomials, while the polynomials 
Vm{x) of order I — 1 or less are known as Van Vleck polynomials. If all were positive, i.e. all Sk and Q negative, 
which can occur for su(l, 1) representations (e.g. boson pairs), then the polynomials Pm{x) would be orthogonal 
polynomials with respect to some weight function. The case of negative residues, rj. < 0, is much less understood and 
a topic of ongoing research^!. 



Appendix C: Continuum limit of the integrable p x + ip y model 

A full analysis of the continuum limit of the integrable p x + ip y has been presented in Rcfs. 22 40, We give here a 
brief overview of the procedure, that applies to the hyperbolic RG model in general, in order to maintain consistency 
with the other sections and to demonstrate that the transitions between the different regions in Fig. [T] arc genuine 
phase transitions and not merely crossovers. The continuum limit of the RG equations (Eq.(17l), is based on the 
following scalings: 

L, M -> 00, (CI) 
M/L = p finite, (C2) 
GL = g finite, (C3) 

Q/L = <l=Yg-\ +P ' (C4) 

|E Sfe /fe) = / e(v)f(v)dn, (C5) 
k A 

~^/(£a) = [ P(x)f(x)dx, (C6) 

where A is the support of the weighted level density g(r]) and f2 the support of the pairon density p(x). Their 
normalization is given by 

' Q(r))dr) = 1, (C7) 

p(x)dx = p. (C8) 
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The RG equations Eq.(17l take the following form in the continuum limit: 



« . 1 r ^)_ dv _ v f k^i^ =0j Vxea (C9) 



x 2 7a x — rj 

The integral over x' is to be understood as a Cauchy principal value, V . To sol ve thi s equation, we start from an 
electrostatic analogy that goes back to work of Stieltjes on orthogonal polynomial d 39 * 41 !, and has been applied in the 
present context by GaudirP^l and was later elaborated in Refs. 26 43 for the rational model and in Refs. I22|4"U1 for the 
hyperbolic model. One considers the pairons to represents M free charges with unit negative charge, positioned in 



the two-dimensional plane at coordinates given by the real and imaginary parts of pairon values E a . Then Eq. ( C9 ) 
represents the two dimensional electrostatic force generated by the repulsion among the free charges, the attraction 
of a fixed charge density g(i]) along the interval A, and the interaction with a central charge q located at the origin, 
derived from a potential V(z) in the two-dimensional complex plane, 

V(z) = -qbx(\z\)-l [ g(rj)ln(\z-r 1 \)dr 1 + [ p(x') ln(\z - x'\)dx'. (CIO) 
2 J a Jn 



Therefore, any particular solution of (C9) gives the equilibrium positions of the pairon density pix) in the two 
dimensional plane. Note that q can change from attractive to repulsive depending on the filling p and on the coupling 
g. This change in sign is associated with a QPT whose interpretation is one of the main goals of this work. To find 




FIG. 12: Contours used to evaluate the continuum limit. 



the unknown pairon density p(x), we assume that it corresponds to a continuous charge distribution along a certain 
curve (or set of curves) tt. The electric field F(z) created by the charge q at the origin and the charge densities p 
(repulsive) and g (attractive) will have a discontinuity across fl given by 

p(x) = F ( x +)- F ( x -) for x e o. (C11 ) 

This allows us to rewrite the Cauchy principal value integral over f2 as an integral along a clockwise contour Cq 
around fi: 

V[**± d x' = ±j ^d, (C12) 



J n x — x' 2ni J Cn x — z 

To model the discontinuity in F(z) we propose the following ansatz: 

F(z) = Rn(z) [ ^d Vl (C13) 

where Rq{z) is a function that explicitly takes into account a cut along fi. If f2 corresponds to a single segment with 
cndpoints a and b, then one can take 

R a - b (z) = y/(z-a)(z-b). (C14) 
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For more segments, Rq(z) has to be a higher order function, e.g. 

n 

Rn(z) = Y[ \/{ z - a i)( z - b i)> 
t=l 

where n is the number of segments. In general, one can write an asymptotic expansion for large z as 

+00 



(C15) 



(C16) 



i=0 



Note that the total charge in the system is given by p — q — 1/2 = —l/(2g). Therefore the asymptotic behavior of F 
is given by 



F(z) = — + O (z- 2 ) , forlzl^oo. 
2gz 



This allows us to put constraints on the unknown function <f)(z) 



4>(i])ri l drj = 0, for i = 0, . . . ,n — 1; 
1 



(C17) 

(C18) 
(C19) 



Now we can rewrite Eq.(C12) by deforming the contour around Q into a contour Ca around A plus a contour C 
around the point at infinity: 



p(x') , , If R n (z) f 4>{r}) 



:dx' 



n x ~ x 
1 



2ni J Cn x - z J A z-T] 



-drjdz 



1 ■ - " ; ' f 7 T7 " \ dzdrj H 

2ttz 7a Jc a ( z ^ - z) 2m 



<l>{n) 



0(77) 



dzdrj 



A a -77 



dry. 



(C20) 



where the integral along Coo vanishes because of Eq.(|C18|) as can be seen in the following expansion: 

Rn 0) 



Coo ( z ~ - z) 



dzdrj 



+<*> n-j-2 



E+oo 
3=0 



3=0 '3' 



+00 



dzdrj 



- V xVj / (j>{rf)r] k dr] </ z n - j - k - l - 2 dz 
j,k,i=o Ja Jc ~ 
+00 „ 
27ri ^ x l rj5k, n -i-j-l / (f>{r])r] k drj 



0. 



(C21) 



because only terms with k < n are allowed by the Kronecker delta. Comparing Eq. ( C20 1 with Eq. ( C9 ) , we find that 
the solution is given by 



q6(rj) + hg(r)) 



The segments of f2 correspond to equipotential curves in the complex plane, determined by 

Re / F(z)dz = 0. 



(C22) 



(C23) 
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The energy density is given by e — E/L and can be evaluated as 

e = {v\Hh\v)/L+ I xp{x)dx 
Jn 

= ~o / w(»7)*7+ — f z Rn(z) / dridz 

2 J A 2m J Ca J\z~r] 

2 



m{ri)dr) + / zR n (z) [ ^-drjdz + / zR n (z) [ ^-drjdz 
2m Jc A Jaz-V 2m /coo J\z-V 

VQ(v)dv + o / VQ(v)dv + 7T~ T2 r j [ 4>{g)v k dr] f 



i=l 



Let us analyze in particular the case of a single segment, i.e. n = 1. Eqs.( |C18|C19| reduce to 



a5 /" ^ ^ 1 1 moo 

dry = -g = - - — - p, (C25) 



2 J A ^( v -a)(r]-b) 2 2g 

T( ^ - (C26) 

It is interesting to compare this to the result of a mean field theory for the ground state (see Appendix [D]) based on a 
grand potential Hh — pS z , with chemical potential p and a pairing field A = G^ fc y^rfkiS^}. One obtains the same 
expressions as in Eqs.( C25|C26 ), provided that one equates 



which means that one can identify 



y/(j)-a){v-b) = Vi.V - 2 V? + 4??A 2 , (C27) 

a + b = 4(p - A 2 ), a& = 4/z 2 (C28) 
a - b = ±4AVA 2 - (C29) 

Again one can distinguish the three regions in the quantum phase diagram: 

A 2 

weak coupling BCS: < =|- < p, and a and 6 complex. The pairons form an arc in the complex plane. 

A 2 

weak pairing at intermediate coupling: < p < and a and b real and negative. The distribution of pairons 
along the negative real axis between a and b has to be supplemented with two arcs in the complex plane that 
close around the origin, touching the real negative axis in o, and a segment of real positive pairons between the 
origin and the point where the arcs touch the positive real axis. 

strong pairing: p < 0, and a and b real and negative. All pairons fall in the segment on the negative real axis 
between a and b. 



For n = 1 one can further elaborate Eq. ( C24 ) for the energy: 



g g 2 L gM V(7 7 -2 At ) 2 + 47 ? A 2 ^ 



where we have defined 



M (2p-l) + ^-6( M ,A 2 ), (C31) 



G(p, a2 ) = \J a g(v)V(v- 2 ^) 2 +^ 2 dv- (C32) 
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We can formally express Eqs.( C25|C26 ) as 

dQ 



elf i 



(p,A')=2p-l, 



d® 
dA 2 



(p,A 2 ) 



(C33) 



Substituting these conditions into Eq.(C31), we see that they guarantee that e(/i, A 2 ) reaches an extremum at the 
corresponding values of p and A 2 , which coincides with the mean field theory of Appendix [d] 

Note that the pairon distribution does not depend on the details of the representation for the su(2) operators S k . 
The only model factor that enters into the equations is the weighted level density q(t]), which corresponds to the 
distribution of Sk in the vacuum state \v). The fact that e, p and g all depend on Q (fi, A 2 ) means that a QPT will 
occur only if Q(p, A 2 ) displays a singularity in one of its derivatives. From Eq.(C32) one can infer that this will be 



the case if the argument of the square root vanishes, which happens at p = if rj = is included in the domain A. 



Appendix D: Mean field theory of the hyperbolic Richardson-Gaudin model 

Mean field theories for separable interactions have been discussed before in the context of p x +ip y pairing models^ 1 ^. 
Here we present a more general derivation applicable to all exactly solvable models derived from the hyperbolic RG 
model, the p x + ip y model representing a particular case. Our mean field theory for the Hamiltonian Hh from Eq.( 16 ) 
is based upon a variational wave function 

\6) =e^ e " s i\v). (Dl) 

For fcrmions, the mean field theory is often presented in terms of a double set of complex parameters Vk/uk = Ok, 
with the condition that \uk\ 2 + \vk\ 2 = 1- These parameters can be related to 9k as 

Ufe = 7rW ^TrW (D2) 



m= Y[( 1 + \e k f) 2s ", S* k \v) = -s k \v}, (D3) 



Let us denote normalized expectation values as (...). From the normalization condition 

k 

one easily finds the normalized expectation values for the basic operators of the algebra: 



I _ Id. |2 



Furthermore, one can show that 

(StSi) = (St) (Si) + 8 k J 8k + W)\ (D5) 



This allows us to write the mean field expectation value for Hh as 



(H h ) = E I* (Sk) ~ G g 0» ( S k) ( S k>) ~ G Y1 1* ^ \f k l))2 - (D6) 

To derive the mean field equations, we will vary the following expression: 

(0\H h \0) - 2p(9\S z \6) - E mf (6\6), (D7) 

where p is a Lagrange multiplier that fixes (S z ) and represents the chemical potential in particle (fermion or boson) 

representations of the model, while E m f is needed to fix the normalization. By varying Eq.( |D7[ ) with respect to (9\ 
we obtain the mean field equation: 

H mf \9)=E mf \6), (D8) 
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with 



(D9) 

(D10) 
(Dll) 



To check whether the ansatz of Eq.(Dl ) is indeed an eigenstate of H m f , we evaluate 



= e E»e*s+ [2^^ + ^^ + ^ [(^A*-A)0h + 2«^] S+-J]V%^ J ( D12 ) 

V fe fc k / 

One sees that |0) is an eigenstate provided that 

(0gA*-AV5fr +20*&=O. 
The corresponding eigenvalue of H m f is given by 

E mf = -2j2 s k [& + 0k y/Vk A*] • 



(D13) 



(D14) 



9k must have the same phase as A, which we can absorb in the definition of the operat or S* , so we can take A and 
Ok to be real without loss of generality. One finds that the ground state solution of Eq.(D13) is given by 



ik = 



?7fcA 



Then, the mean field eigenvalue becomes 

and the expectation values become 

(S z k ) = s k L 



E mf = -2j2skJe k + Vk& 2 , 



(Sk)=s k 



/%A 



tkVk 



A 2 



{Hh) = -^ Sk -yft+^ O 



VkSk 



1 - 



ve k +vkA\ 



Selfconsistency with the definition of A in Eq.( |Dll[ ) gives us the gap equation: 



Vk 



1 



while fixing the expectation value of S z gives us the number equation 



M - L/2. 



In the thermodynamic limit L — > 00, with g — GL finite. The expression for £* simplifies to: 

^ = ri k /2 - /x. 

Adding 1/2 times Eq.(D19l to Eq.( D20[ ) allows us to rewrite the number equation as 



T.Y. 



s k 



= — -1/2 + p. 



(D15) 
(D16) 

(D17) 
(D18) 

(D19) 

(D20) 

(D21) 
(D22) 



2G 



Another interesting situation occurs 



We see that fi = occurs at g = 1 — 2p, i.e. special case (i) from Sec. IIB 
when /i = Tj-, because then one has that \Jt%. + f?fcA 2 = \r)k/2 + fj,\ and consequently 6k = A/fc. By adding 1/2 times 
Eq.(D19) to Eq.(D22) for this case, one obtains (assuming all rjk > 0) that g~ x = 1 — p, i.e. special case (ii) from 
Sec.LLLB 
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